Antiferromagnetic order in systems with doublet Stot = 1/2 ground states 



O 

00 



0-) 

i: 



o 



> 

O 



X: 



Sambuddha Sanyal,^ Argha Banerjce/ Kedar Damle/ and Anders W. Sandvik^ 

^Department of Theoretical Physics, Tata Institute of Fundamental Research, Mumbai 400005, India. 
^Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachussetts 02215, USA. 

We use projector Quantum Monte-Carlo methods to study the Stot = 1/2 doublet ground states 
of two dimensional S = 1/2 antiferromagnets on a L x L square lattice with an odd number of sites 
A'^tot = ■ We compute the ground state spin texture $^(r) = (S^(r))-f in |G)-t-, the Sfot = 1/2 
component of this doublet, and investigate the relationship between , the thermodynamic limit 
of the staggered component of this ground state spin texture, and m, the thermodynamic limit of 
the magnitude of the staggered magnetization vector of the same system in the singlet ground state 
that obtains for even A'tot. and m would have been equal if the non-zero value of Stot in 1^)-^ 
caused the direction of the staggered magnetization vector to be fully pinned in the thermodynamic 
limit. By studying several different deformations of the square lattice Heisenberg antiferromagnet, 
we establish that this is not the case. For the sizeable range of m accessed in our numerics, we 
find a univeral relationship between the two, that is independent of the microscopic details of the 
lattice level Hamiltonian and can be well approximated by a polynomial interpolation formula: 
"-^ ~ (|-f-i)™- + + ''"i^' " ~ 0-288 and b ^ -0.306. We also find that the full 

spin texture <l?^(r) is itself dominated by Fourier modes near the antiferromagnetic wavevector 
in a universal way. On the analytical side, we explore this question using spin-wave theory, a 
simple mean field model written in terms of the total spin of each sublattice, and a rotor model for 
the dynamics of n. We find that spin- wave theory reproduces this universality of "l?^(r) and gives 
= {l-a- p / S)m+{a/ S)m^ + 0{S''^) with a ^ 0.013 and /3 ^ 1.003 for spin-S antiferromagnets, 
while the sublattice-spin mean field theory and the rotor model both give = for S = 1/2 
antiferromagnets. We argue that this latter relationship becomes asymptotically exact in the limit 
of infinitely long-range unfrustrated exchange interactions. 



PACS numbers: 75.10.Jm 05.30.Jp 71.27.+a 



I. INTRODUCTION 

Computational studies of strongly correlated systems 
necessarily involve an extrapolation to the thermody- 
namic limit from a sequence of finite sizes at which cal- 
culations are feasible. Understanding^ and at times 
reducing^ these finite-size corrections to the thermody- 
namic limit is thus an important aspect of any such cal- 
culation. For instance, the best estimates of m, the mag- 
nitude of the ground state Neel order parameter in the 
thermodynamic limit of the two-dimensional 5 = 1/2 
square lattice Heisenberg antiferromagnet rely on a se- 
quence of Lx X Ly systems with even length (Ly) in the 
X (y) direction and periodic boundary conditions in both 
directionsi^i^ Other studies suggest^ that it is some times 
advantageous to use "cylindrical" samples with periodic 
boundary conditions in one direction and pinned bound- 
ary conditions in the other direction, whereby spins are 
held fixed by the use of pinning fields on one pair of 
edges — this choice also allows for a very accurate deter- 
mination of ground-state parameters such as m for spe- 
cific values^ of the aspect ratio Ly/L^. 

All these approaches focus on systems with an even 
number of spin-half variables; this choice allows the 
ground-state of the finite system to lie in the sin- 
glet sector favoured by unfrustrated antiferromagnetic 
interactions 1^ Although not commonly used, another 
choice is certainly possible: Namely, one could in princi- 
ple consider antiferromagnets on a L x L square lattice 
with an odd number A^tot = L" of spin-half moments. 



Such a system is expected to have a doublet ground state 
with total spin Stot = 1/2. Focusing on the S^^^. = 1/2 
member |G)^ of this doublet, one could examine the 
ground state spin texture defined by = (55)^ 

(where (...)-[- refers to expectation values in \G)-f), and 
use the antiferromagnetic component of this spin texture, 
defined as 
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At, 



(1) 



to obtain information about the antiferromagnetic order- 
ing in the system (here = +1 on the A sublattice and 
— 1 on the B sublattice). 

Clearly, provides a measure of antiferromagnetic 
order that is quite distinct from the conventional order 
parameter m, which can be defined, e.g., according to 
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where (• . • )o denotes averages in the singlet ground state 
realized for even A'tot- The relationship between the 
thermodynamic-limit values of and m is a fundamen- 
tal aspect of the spontaneously broken SU{2) symmetry 
of the Neel state. However, not much is known about 
it beyond the fact that is significantly smaller than 
m for the nearest neighbour Heisenberg antiferromagnet 
on ths square latticei^ Here, we provide a more detailed 
characterization of this relationship. 
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Our basic result is that is determined in a universal 
way by the value of m. In other words, plotted against 
m for several different deformations of the S* = 1 /2 square 
lattice Heisenberg antiferromagnet falls on a single curve 
which defines a universal function that is insensitive to 
the microscopic details of the model Hamiltonian. This 
universal function is well- approximated by a polynomial 
interpolation formula: 

n"^ « (i - ^ - + am^ + 6m^, (3) 

with a « 0.288 and b w —0.306. In addition, we also find 
that the full spin texture $^(r) is dominated by Fourier 
modes near the antiferromagnctic wave-vector in a uni- 
versal way independent of microscopic details. We show 
that this universality is captured by spin-wave theory, 
which also predicts 

= {l-a- I3/S)m + {a/S)m'^ + 0(5'-2), (4) 

with a « 0.013 and /3 « 1.003 for spin-S* antifcrromag- 
nets. In addition, we explore two other ways of thinking 
about this universal function. One of them is a mean 
field theory formulated in terms of the total spin of each 
sublatticc, while the other approach is in terms of a quan- 
tum rotor Hamiltonian for the Neel vector n of a system 
with an odd number of sites. Both these give 



for S = 1/2 antiferromagnets, which is close to the 
observed relationship but not exactly right. We argue 
that this latter estimate (Eqn. [S|) will become asymp- 
totically exact in the limit of infinitely long-range un- 
frustrated exchange interactions. In this limit, we also 
expect TO — > 1/2, and our polynomial fit to the universal 
function n^(m) was therefore constrained to ensure that 
—?' m/3 when to — >■ 1/2. 

The outline of the rest of the paper is as follows: In Sec- 
tion |ll] we define various deformations of the square lat- 
tice 5 = 1/2 Heisenberg antiferromagnet. In Section Hill 
we outline the projector quantum Monte Carlo (QMC) 
method used in this study, and then discuss in some de- 
tail our QMC results for as well as the full spin tex- 
ture ^^(r), focusing on the universal properties alluded 
to earlier. In Section ITVl we outline three analytical ap- 
proaches to the relationship between and to. The first 
is a large-S* spinwave expansion, within which we calcu- 
late the ground state spin texture $^(r) and its antiferro- 
magnctic Fourier component to leading 0{l/S) order, 
and demonstrate that such a calculation also yields the 
universality properties summarized earlier, but does not 
provide a quantitatively accurate account of the QMC 
results for <I>^(r) or n^{m). The second is a mean-field 
theory formulated in terms of the total spin of each sub- 
lattice. And the third approach is in terms of a quantum 
rotor Hamiltonian which is expected to correctly describe 
the low-energy tower of states for odd A^tot- In ScctionlVl 
we conclude with some speculations about a possible ef- 
fective field theory approach to the calculation of <I>^(f'). 
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FIG. 1: An illustration of the interactions present in the J J' 
(left panel) and JJ2 (right panel) model Hamiltonians. In this 
illustration, black bonds denote exchange interaction strength 
of J, while a red bond represents exchange strength of J' ( J2) 
in the left (right) panel 
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FIG. 2: Bond and plaquette operators in JQ model Hamil- 
tonians. A thick bond denotes a bipartrite projector acting 
on that bond. All possible orientations of these bond and 
plaquette operators are allowed. 

II. MODELS 

We consider four deformations of the square lattice 
5=1/2 nearest neighbour Heisenberg antiferromagnet; 
all four retain the full SU{2) spin rotation symmetry of 
the original model. 

The first of these models is the couplcd-dimer antifer- 
romagnet, in which there are two kinds of nearest neigh- 
bour interactions J and J', as shown in Fig. [1] (left panel), 
where the ratio a ~ J' /J can be tuned from a = 1 to 
a = ttc ~ 1.90 at which coUinear antiferromagnctic order 
is lostii The Hamiltonian for this system reads: 

Hj.r = J ^ S, • + J' ^ S, • S, , (6) 

<y> (ij)' 

where (ij) {{ij) ) denotes a pair of nearest neighbour sites 
connected by a black (red) bond (see Fig. [1}. Another 
deformation of the Heisenberg model, the JJ2 model, has 
additional next nearest neighbour Heisenberg exchange 
interactions J2, as shown in Fig. [1] (right panel). The 
Hamiltonian reads 

Hjj, = S, • S, + J2 ^ S, • S„ (7) 

(ij) ((y)) 

where {{ij}) denotes a pair of next nearest neighbour 
sites. Both these are amenable to straightforward spin- 
wave theory analyses, and the coupled dimer model can 
also be studied numerically to obtain numerically exact 
results even for very large sizes due to the absence of any 



sign problems in Quantum Monte Carlo studies. How- 
ever, exact numerical results on the JJ2 model are re- 
stricted to small sizes since Quantum Monte Carlo meth- 
ods encounter a sign problem when dealing with next- 
nearest neighbour interactions on the square lattice. 

In addition, we study two generalizations that involve 
additional multispin interactions; the " JQ" modelsjSii^ Of 
these, the JQ2 model has 4-spin interactions in addition 
to the usual Heisenberg exchange terms, and is defined 
by the Hamiltonian 
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JQ2 



E 



(8) 



where the plaquettc interaction Q2 involves two adjacent 
parallel bonds on the square lattice as shown in Fig. [2] 
(middle panel) and 



p — 1 _ s • S 



(9) 



is a bipartite singlet projector. The first term in Eqn. [8] 
is just the standard Heisenberg exchange. Similarly, the 
JQ3 model has 6-spin interactions and is defined by the 
Hamiltonian 



H 



JQ3 
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(ij^kl^nm) 
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(10) 



where the plaquette interactions now involve three ad- 
jacent parallel bonds on the square lattice, as shown in 
Fig. [2] (right panel). The products of singlet projectors 
making up the Q2 and terms tend to reduce the Neel 
order of the ground state, and, when sufficiently strong, 
lead to a quantum phase transition into a valence-bond- 
solid statcj^i^ Here we stay within the Neel state in both 
models, and study universal aspects of this state as the 
Neel order is weakened. 



III. PROJECTOR QMC STUDIES 

We use the total spin-half sector versioni^ of the 
valence-bond basis projector QMC methodiiii^ to study 
L X L samples with L odd and free boundary conditions. 
We compute ^^{r) and in such samples for the J J' 
model and JQ models in their antiferromagnetic phase. 
We also study the same models on L x L lattices with L 
even and periodic bondary conditions using the original 
singlet sector valence bond projector QMC method. In 
both cases we use the most recent formulation with very 
efficient loop updatesj ^°i^^ Our system sizes range from 
L = 11 to L = 101, and projection power scales as 
to ensure convergence to the ground state. We perform 
> 10^ equilibration steps followed by > 10^ Monte Carlo 
measurements to ensure that statistical and systematic 
errors are small. 

Data for from a sequence oi L x L systems with L 
odd shows that extrapolates to a finite value in the 
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FIG. 3: An illustrative example of finite size corrections of 
and m^, observed in the antiferromagnetic phase of the J J' 
model(J' — 1.8). Note the non-monotonic behaviour of finite 
size corrections for , which is fitted to a cubic polynomial. 
In contrast, finite size data for is well described by a linear 
dependence on 1/L. 



L — > 00 limit as long as the system is in the antiferro- 
magnetic phase. However, we find that the approach of 
this observable to the thermodynamic limit has a non- 
monotonic behaviour. To obtain accurate extrapolations 
to infinite size, it is therefore necessary to fit the finite size 
data to a third-order polynomial in 1/L. Wc find that the 
coefficient for the leading 1 /L term in this polynomial is 
rather small; this is true for all the models studied here, 
as long as they remain in the antiferromagnetic phase. In 
Fig.. [3] and Fig.[31 we show examples of this behaviour of 
the finite size corrections in . In these figures, we also 
show the approach to the thermodynamic limit for m, as 
measured in a sequence of periodic L x L systems with L 
even. We find that in complete contrast to the behaviour 
of n^, TO extrapoloates monotonically to the thermody- 
namic limit, with a dominant 1/L dependence — this is 
consistent with previous studies of the structure factor 
in square lattice antiferromagnetsi^ (however, with spa- 
tially anisotropic couplings, one can also observe strong 
non-monotonicity in m^). 

The non-zero value of 71^ in the thermodynamic limit 
clearly reflects the long-range antiferromagnetic order 
present in the system and a partial breaking of the SU(2) 
symmetry (due to the fact that we study only one mem- 
ber of the doublet ground state). For periodic systems, 
the same long range antiferromagnetic order is captured 
by the non-zero value of to in the large L limit — and a 
calculation of to (through (m^)) for the odd-L systems 
with periodic boundaries would of course lead to the same 
value. However, since to 7^ n^, the full staggered mag- 
netization is not forced to lie along the z spin axis, and 
it is interesting to ask: What is the relationship between 
these two measures of antiferromagnetic order? Our nu- 
merical data are unequivocal as far as this relationship is 
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a = 0.288 ± 0.022, b = -0.306 ± 0.03 1 
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FIG. 4: Another illustrative example of finite size corrections 
of and m^, observed in the antiferromagnetic phase of JQ2 
model at Q2 = 1.0. Again, note the non-monotonic behaviour 
of finite size corrections for n^, which is fitted to a cubic 
polynomial (only L > 20 data used in the fit). In contrast, 
finite size data for is well described by a linear dependence 
on 1/L. 



concerned, as is clear from Fig. [SJ which shows a plot of 
versus m in the thermodynamic limit of the JJ', JQ2 
and JQ3 models. Here each point represents the result 
of a careful extrapolation similar to the examples shown 
in Fig. [3] and Fig.|4l and provides an accurate estimate of 
the corresponding thermodynamic limits for and m. 
From this figure, it is clear that is a universal func- 
tion of m independent of the microscopic structure of the 
Hamiltonian. To model this universal function, we use a 
polynomial fit that is constrained to ensure that — )■ y 
when TO — i; the rationale for this constraint will be- 
come clear in Sec. |lVl We find (Fig. [5|) that the QMC 
results for n^[m) are fit well by the following functional 
form: 

{m) ^ - ^ - + arri^ + bm^ , (11) 

with a « 0.288 and b « -0.306. 

If one views this universal relationship as being a prop- 
erty of the low energy effective field theory of the anti- 
ferromagnetic phase, one is led to expect that the full 
spatial structure of the spin texture (r ) should also be 
universal. More precisely, one is led to expect that this 
texture is dominated in a universal way by Fourier com- 
ponents near the antiferromagnetic wavevector. To test 
this, we compare the spin texture in the J J' model and 
the JQs model, choosing the strengths of the J' interac- 
tion and the Q3 interaction so that both have the same 
value of TO, and therefore the same value of n^. This is 
shown in Fig. |6l which shows that these very different mi- 
croscopic Hamiltonians have spin-textures whose Fourier 
transform falls on top of each other at and around the 
antiferromagnetic wavevector. 



FIG. 5: Extrapolated thermodynamic values of for three 
different models of antiferromagnets on an open lattice, plot- 
ted as function of staggered magnetisation m for the same 
models on periodic lattices. The former is clearly an uni- 
versal function of the later. This universal function can be 
well approximated by a polynomial fit constrained to en- 
sure that n^(m) — >■ m/3 in the limit of m — ^ |: ~ 
(1/3 - a/2 - 6/4)m + am^ + bm^ , with a ^ 0.288 and 
b ^ -0.306. 



IV. ANALYTICAL APPROXIMATIONS 

We now present three distinct analytical approaches 
to understanding these numerical results presented in the 
previous section: First, we develop a spin- wave expansion 
that becomes asymptotically exact for large S^. Sec- 
ond, we explore a mean-field theory written in terms of 
the total spin of each sub-lattice. Finally, we describe an 
alternative approach in which the low-energy antiferro- 
magnetic tower of states of a spin- 1/2 antiferromagnet is 
described by a phenomenological rotor model^^ adapted 
to the case of a system with odd iVtot- 



A. Spin- wave expansion 

The leading order spin-wave calculation proceeds as 
usual by using an approximate representation of spin 
operators in terms of Holstein-Primakoff bosons. The 
resulting bosonic Hamiltonian is truncated to leading 
(quadratic) order in boson operators to obtain the first 
quantum corrections to the classical energy of the system. 

As is standard in the spin wave theory of Neel ordered 
states, we start with the classical Neel ordered config- 
uration with the Neel vector pointing along the z axis, 
which corresponds to Sp = rjpS. We then represent the 
spin operators at a site r of the square lattice in terms 
of canonical bosons to leading order in S as follows: For 
sites r belonging to the A sublattice we write 

Sj = V2Sbp- SZ = S~blbp, (12) 
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while on sites r belonging to the B sublattice wc write 
SZ = V2Sbp ; Sp ^-S + bibp . (13) 

The number of bosons at each site thus represents the 
effect of quantum fluctuations away from the classical 
Neel ordered configuration. 

To quadratic order in the boson operators, this ex- 
pansion yields the following spin wave Hamiltonian in 
the general case (with arbitrary two-spin exchange cou- 
plings) : 

Hsyj = eciS^ + ^hH'Ih , with 

^' = (4:; J:;) 
w=(;;r). (14) 

Here edS^ is the classical energy of the Neel state, M 
in the first line is a 2A'tot dimensional matrix specified 
in terms of A'tot dimensional blocks A and B, and b is 
a 2A'tot dimensional column vector as indicated above. 
Elements of A and B can be written explicitly as 

App, = {Z^I ^Zf)Spp+Jfp,, (15) 
Bpp, = Jllp. (16) 

In the above, jSp, are Hcisenberg exchange couplings be- 
tween two sites r and r' belonging to the same sub-lattice, 
Jpp, are the Heisenberg exchange couplings between sites 
belonging to different sublattices, and 

Zp = ^J^,, (17) 

r' 

Zp = '^jfp,. (18) 
f' 

The effects of quantum fluctuations on the classical 
Neel state can now be calculated by diagonalizing this 
Hamiltonian by a canonical Bogoliubov transformation 
S which relates the Holstein-Primakoff bosons b to the 
bosonic operators 7 corresponding to spin-wave eigen- 
states 

b = §r, ^, = (j,l) > (19) 

where § is a 2A'tot dimensional matrix that transforms 
from b which creates and destroys bosons at specific lat- 
tice sites r to r which creates and destroys spin- wave 
quanta in specific spin-wave modes ^. Naturally, we must 
require that iJ^ui be diagonal in this new basis. We rep- 
resent this diagonal form as 

= e,iS^- + f rtOT , (20) 

where 



with A denoting the diagonal matrix with the TVtot posi- 
tive spin wave frequencies on its diagonal. 

To construct a § that diagonalizes Hgw in the T basis, 
we look for 2A'^(:,ot dimensional column vectors 

- ' (22) 

which satisfy the equation 

My^ = e^Jy^ (23) 

with positive values of equal to the positive spin-wave 
frequencies for n = 1,2,3... A^tot- Here and are 
A^tot dimensional vectors, 

and 1 is the A'tot x A'tot identity matrix. With these 
in hand, one may obtain A^tot additional solutions 
to Eqn. I23i this time with negative eNt^t+t^ ~ 
interchanging the roles of the A^tot dimensional vectors 
and Up in this construction. In other words, we have 

?/^-+^ - {S) ' (25) 

with ^ = 1,2,3... A^tot- 

We now construct S by using these (with /i = 
1,2,3... 2A^tot) as its 2A'tot columns: 

(y\y2,y3 ..y^^'tot) (25) 

Clearly, this choice of § satisfies the equation 

MS = ISID (27) 

Furthermore, the requirement that the Bogoliubov trans- 
formed operators 7 obey the same canonical bosonic com- 
mutation relations as the b operators implies that § must 
satisfy 

Si'TS = I , (28) 

This constraint is equivalent to "symplectic" orthonor- 
malization conditions: 

iu^)^u-' ~ {v'')^,'' = 5,, , (29) 
(uM)t„^ - (i;'^)tu'' = , 

for II, V = 1, 2, 3 . . . A^tot- It is now easy to see that EgnlTTl 
and Eqn [28] guarantee that is indeed diagonal in the 
new basis, since 

b^Mb = r^§t7\,/§r = r^s^isiDT = rt£ir . (30) 

For periodic samples, it is possible to exploit the trans- 
lational invariance of the problem and work in Fourier 
space to obtain these spin-wave modes and their wave- 
functions and calculate to = S' — A correct to leading 




FIG. 6: Fourier transform (with antiperiodic boundary conditions assumed for convenience) of the numericaUy computed (for 
J J' and JQs model with L = 65, S = 1/2 ) [r) along cuts passing through the antiferromagnetic wavevector (tTjTt). Note 
the universality of the results in the neighbourhood of the antiferromagnetic wavevector, which in any case accounts for most 
of the weight in Fourier space. 



order in the spin-wave expansion — as these results are 
standard and wcll-knownfi^ we do not provide further 
details here. On the other hand, the corresponding re- 
sults for L X L samples with free boundary conditions 
and Na = Nb -I- 1 do not seem to be available in the 
literature, and our discussion below focuses on this case. 

We begin by noting that the non-zero entries in A only 
connect two sites belonging to the same sublattice, while 
those in B always connect sites belonging to opposite 
sublattices. As a result of this, the solutions to the equa- 
tion for can also be expressed in terms of a single 
function /^(r) defined on sites of the lattice. To see this, 
we consider an auxiliary problem of finding such that 
the operator A — B ~ e^r]f has a zero mode f^{r) (as 
before, ry? is -1-1 for sites belonging to the A sublattice, 
and —1 for sites belonging to the B sublattice). 

This auxiliary problem has A'tot solutions correspond- 
ing to the A'tot roots e^i of the polynomial equation 
det(A — B — e^?7r) = 0; these can be of cither sign. To 
make the correspondence with the positive e^^ solutions 
(«'',«'') (with ^ = l,2...A'tot) of the original equation 
My^ = e^Iy^, we now note that 

{f,,\A-B\f^) = ~e,,N^ (31) 

where 

iV^^^|/,,(r^)p-^|/^(rB)p. (32) 

rA TB 

Since A — B is a positive (but not positive definite) oper- 
ator, this implies that has the same sign as N^^ for all 
non-zero e^. To make the correspondence with the pos- 
itive = solutions {fi = 1,2... iVtot) of the original 



problem, we can therefore make the ansatz 

=/,.(rA)/v^,< =0 (33) 
< =-/m(^b)/\/^,<, =0 

if Nfj, > 0, or the alternative ansatz 

Ks^-U^b)/ <a =0 (34) 

< =/M(^A)/y^,< =0 

if Nf^ < 0. Here, ta (tb) denotes sites belonging to 
the A (B) sublattice of the square lattice. This ansatz 
clearly ensures that the (with /i = 1, 2..iVtot) obtained 
in this manner satisfy the original equation with positive 

= and are appropriately normalized. 

Atlhough this approach is not the one we use in our ac- 
tual computations (see below) , it provides a useful frame- 
work within which we may discuss possible zero frequency 
spin- wave modes, i.e X^g = for some A mode 
with A^,-| = clearly corresponds to a putative zero eigen- 
value of the operator A — B. From the specific form of 
A — B in our problem, it is clear that such a zero eigen- 
value does indeed exist, and f^ai^^ the corresponding 
eigenvector oi A — B, can be written down explicitly as 

/mo(^1 = 1 (35) 

Since this corresponds to the root e^^ = of the auxiliary 
problem, it can in principle be used to obtain a pair of 
zero frequency modes e^^ and e^o-i-A'tot foi' the original 
problemof finding and y^ that satisfy My'^ = e^Iy^ . 

However, we need to ensure that the symplectic or- 
thonormalization conditions (Eqn. I30p arc satisfied by 
our construction of the corresponding y^^^' and yAio-i-JVtot ^ 
This is where the restriction to a A^tot ~ L x L lattice 
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with Na = Nb + 1 enters our discussion. For this case, 
N^g = Na — Nb = 1, and we are thus in a position 
to write down properly normahzed zero-mode wavefunc- 
tions: 



where 



Mo 



= /mo(^a),< =0 



and 



yJVtot+MO ^ yM 



(36) 



(37) 



[ParentheticaUy, we note that the question of zero fre- 
quency spinwave modes for the more famihar case with 
A^^ = Nb and periodic boundary conditions has been 
discussed earher in the hteraturoi^ and will not be con- 
sidered here.] 

Thus, the equation Aly^ = e^jly^ has a pair of zero 
modes related to each other by interchange of the u and 
V components of the mode, and it becomes necessary to 
regulate intermediate steps of the calculation with a stag- 
gered magnetic field zehTjp with infinitesimal magnitude 
e/i > in the z direction. Denoting the corresponding A 
by yl'^'' , we see that A^*^ — B is now a positive definite 
operator and does not have a zero eigenvalue. Indeed, 
it is easy to see from the foregoing that the correspond- 
ing eigenvalue now becomes non-zero, yielding a positive 
spin-wave frequency AJ^'j^ = A^totf/i- One can also calcu- 
late the C(e/i) term of /^J;(?^) and check that /^[j tends 
to fiigir) in a non-singular way as Ch — 0, from which 
one can obtain the corresponding y^°{£h) analytically in 
this limit. Thus, the contribution of the zero mode to all 
physical quantities can be obtained in the presence of a 
small e/i > 0, and the — > limit of this contribution 
can then be taken smoothly and analytically at the end 
of the calculation. 

In our actual calculations, we use this analytical un- 
derstanding of the zero frequency spin wave mode to 
analytically obtain the properly regularized zero mode 
contribution to various physical quantities, while using a 
computationally convenient approach to numerically cal- 
culate the contribution of the non-zero spin wave modes. 
To do this, we rewrite Eqn. [23]for ^ = 1,2,3... A'tot as 



where 



{A + B)r = 

{A-B)V = 



(38) 



(39) 



1P^' 

This implies 

{A-B){A + B)^^ 
{A + B){A- B)yj^' = Xf,iA + B)^' 
We now decompose 

A-B = K^K. 



X,{A-B)r = Xlr (40) 



(41) 



(42) 



(43) 



with uj the diagonal matrix with diagonal entries given 
by eigenvalues of the real symmetric matrix A—B, and U 
the matrix whose rows are made up of the corresponding 
eigenvectors. 

With this decomposition, we multiply Eqn ST] by K 
from the left to obtain 



K{A + B)Kh''^Xf,x^ 



(44) 



with = Ktp^. From the solution to this equation, we 
may obtain the cf) as 



(45) 



and thence obtain ip^ using Eqn 1391 In order to ensure 
the correct normalization of the resulting u^^v^, we im- 
pose the normalization condition 



(46) 



Thus our computational strategy consists of obtaining 
eigenvalues of the symmetric operator K{A + B)K\ and 
using this information to calculate the and thence 
the Bogoliubov transform matrix §. Notwithstanding the 
normalization used in Eqn 1461 the zero mode with = 
causes no difficulties in this approach, since we work in 
practice with the projection of K{A + B)K'^ in the space 
orthogonal to the zero mode. This is possible because 
we already have an analytic expression correct to 0{eh) 
for y^°{eh) and y^^°^^^° {^h) corresponding to this zero 
mode, and do not need to determine these two columns 
of S by this computational method. 

We use this procedure to calculate the zero tempera- 
ture boson density as 



(47) 



In this expression, one may use the numerical procedure 
outlined above to obtain the contribution of all ^ 
directly at = 0, while being careful to use our analyti- 
cal results for v^°{eh) to obtain the limiting value of the 
contribution from /.( = yUo- This gives 



M#MO 



(48) 



(49) 



Here, the distinction between sites on the A and 
B sublattices arises in this final result because 
-1 for r belonging to the B sublat- 



lim^^^o v't°{eh) 



tice, while Wui^^^q Vp° [eh] ~ for belonging to the A 



sublattice. 



Knowing the average boson number at each site gives 
us the first quantum corrections to the ground state ex- 
pectation value {S^{r)): 



ag= 0.476, bg= 0.50000, Cg= -0.00038, dg= 0.85 



(5^(f))=7y,.(5-(6t6,~)) 



(50) 



This result for the spin-wave corrections to the ground 
state spin texture then allows us to write n 



= S-A 



(51) 



where A represents the leading spin-wave correction to 
the classical value for n^. 

In order to obtain reliably in this manner, it is im- 
portant to understand the finite size scaling properties 
of A for various values of J / J in the striped interaction 
model and J2/J in the model with next-nearest neigh- 
bour interactions. In Fig. [71 wc show a typical exam- 
ple of this size dependence. As is clear, we find that A 
has a non monotonic dependence on L: A initially in- 
creases rapidly with size, and, after a certain crossover 
size L*, it starts decreasing slowly to finally saturate to 
its asymptotic value. This non-monotonic behaviour is 
qualitatively similar to that observed in the finite size 
extrapolations of from our QMC data earlier. To ex- 
plore this unusual size dependence further and reliably 
extrapolate to the thermodynamic limit, we analyze the 
contributions to A from the spin-wave spectrum in the 
following way: We note that there is always a mono- 
tonically and rapidly convergent 0(1) contribution to A 
from the lowest frequency spin-wave mode, whose spin- 
wave frequency scales to zero as l/A'tot (for any finite 
iVtot, this is not an exact zero mode of the system). We 
dub this the 'delta-function contribution' and its thermo- 
dynamic limit is easy to reliably extrapolate to. In ad- 
dition, there is a 'continuum contribution' coming from 
all the other spin- wave modes, each of which contributes 
an amount of order ©(l/iVtot). This contribution con- 
verges less rapidly to the thermodynamic limit, and also 
happens to be non-monotonic: it first increases quickly 
with increasing size, and then starts decreasing slowly to 
finally saturate to the thermodynamic limit. 

The delta-function contribution can be fit best to a 
functional form 



FsiL) 



, C5 as 



Z3' 



(52) 



with the dominant 1/i^ term accounting for the mono- 
tonic increase with L, while the continuum contribution 
is fit to 



Fc{L) 



Cc 

L 



ac_ 
L2 



L3' 



(53) 



whereby the size dependence is predominantly deter- 
mined by the competition between the term proportional 
to \/L which decreases with increasing L, and the term 
proportional to l/L^ which increases with increasing L. 
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FIG. 7: A typical example of the finite size scaling of the delta- 
function and continuum contributions to A. Note the mono- 
tonically increasing size dependence of the delta- function con- 
tribution, and the non-monotonic and more slowly converging 
nature of the continuum contribution. Due to this difference 
in their behaviour, we find it more accurate to separately fit 
each of these contributions to a polynomial in 1/L and use 
this to obtain the thermodynamic limit of the total A. Here 
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This gives rise to non-monotonic behaviour whereby the 
continuum contribution first increases rapidly and then 
decreases slowly beyond a crossover length L* to finally 
saturate to its infinite volume limit. We also find that the 
length L* gets larger as we deform away from the pure 
square lattice antiferromagnet, making it harder to ob- 
tain reliable extrapolations to the thermodynamic limit. 

Using such careful finite-size extrapolations to obtain 
A for various values of J2/J and J / J, we compare the 
result with A' calculated analytically. Specifically, we 
now ask if the universality seen in our QMC results is re- 
flected in these semiclassical spin-wave corrections to 
and m. The answer is provided by Fig. |8l which shows 
that the numerically obtained spin-wave corrections ap- 
parently satisfy a universal linear relationship 



A - A' ss 1.003 + 0.013A' 



(54) 



as one deforms away from the pure square lattice anti- 
ferromagnet in various ways. 

What does this imply for n^im) to leading order in 
1/5? To answer this, we note that 



(55) 



Using our numerically established universal result to re- 
late A — A' to A' and thence to m itself, we obtain the 
universal relationship 



(56) 



with a w 0.987 - imZ/S and P « COIB/S*. However, 
being a large-S* expansion, spin-wave theory is unable to 
give a quantitatively correct prediction for (to) for the 
S — 1/2 case. 

Finally, we use our spin-wave predictions for the 
ground-state spin texture to look at the Fourier trans- 
form of the spin-texture for various deformations of the 
pure antiferromagnet. The results arc shown in Fig. [HI 
which demonstrates that spin-wave theory also predicts 
that the Fourier transform of the spin-texture near the 
antiferromagnetic wave-vector is a universal function of 
the wavevector; this provides some rationalization for the 
observed universality of the Fourier transformed spin tex- 
ture seen in our QMC numerics. 



quantum number of Sa is Sb + 1/2 while the total spin 
quantum number of Sb should be taken to be Sb, where 
Sn ~ No 10. tends tn infinitv in the thermndvnamic limit. 
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A-A'= 1.003 + 0.013 A' 
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FIG. 8: A — A', the difference between the leading spin wave 
corrections to and m, plotted against the leading spin- wave 
corrections A' to m for the J J and JJ2 models described in 
the text. 



Hamiltonian that describes the low energy part of the 
spectrum: 



MF 



— JufSa ■ Si 



(57) 



with Jmf > 0. Within this mean-fleld treatment, the 
5tot = 1/2, 5*^04 = 1/2 ground state that we focus 
on in our numerics is thus the ^tot = 1/2, 5*^ = 1/2 
state obtained by the quantum mechanical addition of 
angular momenta Sb and Sb + 1/2. Within this mean- 
field theory, is modeled as the expectation value of 
(5^ — Sg)/Ntot in this state, which can be readily ob- 
tained in closed form using the following standard re- 
sult for the minimum angular momentum state |J ~ 
ji — j2,TOj) state obtained by the addition of angular 
momenta ji and 72 (with ji > J2): 



(ii , mi ; j2 , ^2 1 J, TO j) = pjc 



(58) 



B. Sublattice-spin mean-field theory 

We now turn to a simple mean-field picture in terms 
of the dynamics of the total spins Sa and Sb of the A 
and B sublattices respectively. When Na = Nb + 1, 
it is clearly appropriate to assume that the total spin 



with 



and 



PJ 



'(2J+l)!(2j2)! 



(2ji 



D! 



(59) 



J 



[(^-^ + TOi)!((ji - TOi)!]^/' [(j2 + m2)!(j2 - m^W + mj)\{J - m,,)!]"^/' 



(60) 
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for mi + 7712 — mj and c;^'-,^:^^ = otherwise. 

In our case, ji ^ Sb + 1/2, j2 = Sb, J = 1/2, 
TO,/ = 1/2, and = (mi — m2)j,mj/^tot can therefore 
be readily calculated to obtain 



1^3 + l) /iVtot 



(61) 



within this phenomenological approach. 

On the other hand, when iV^ = Nb, we may also 



calculate m- 



{{Sa — Sb) )j=o/N^ot within the same 



sublattice-spin approach 
m^ (451 



-^Sb)/N, 



tot 



(62) 



This allows us to compute the ratio /m in the thermo- 
dynamic limit: 



Is there a limit in which this sublattice-spin mean-field 
theory is expected to give exact results? To answer this, 
we note that the sublattice-spin model represents the 
Hamiltonian of an infinite-range model in which every A 
sublattice-spin interacts with every B sublattice-spin via 
a constant (independent of distance) antiferromagnetic 
exchange coupling Jmf- Thus, our mean- field theory is 
expected to become asymptotically exact in the limit of 
infinitely long-range unfrustrated couplings. In this limit, 
we also expect m — 1/2, and thus, our mean field theory 
predicts that m/3 when m 1/2. This is the 

constraint that we built into our choice of polynomial fit 
for n^{m) in Sec. IIIIl 



C. Quantum rotor Hamiltonian 

When any continuous symmetry is broken, the corre- 
sponding order parameter variable becomes very "heavy" 
in a well-defined senseJ^ The long-time, slow dynamics 
of this heavy nearly classical variable is controlled by 
an effective "mass" that diverges in the thermodynamic 
limit. 

For a Neel ordered magnet, the order parameter is the 
Neel vector n. In the usual case of an antiferromagnet 
with an even number of 5 = 1/2 moments, the low- 
energy effective Hamiltonian that controls the oricnta- 
tional dynamics of the Neel vector n is 



L- L 

^xNtot 



(64) 



where L is the angular momentum conjugate to the 
"quantum rotor" coordinate h = n/\n\, x is the uniform 
susceptibility per spin, and A^tot is the total number of 
spins. 

What about our case with Na = Nb + 1 and an odd 
number of spins A'tot? Following earlier work on quantum 



rotor descriptions of insulating antiferromagnets doped 
with a single mobile charge-carrier—, we postulate that 
the correct rotor description of our problem is in terms of 
a rotor Hamiltonian in which L is replaced by the angular 
momentum operator L' conjugate to a quantum rotor 
coordinate n that now parametrizes a unit-sphere with a 
fundamental magnetic monopole at its originJ^ In other 
words, we postulate a low-energy effective Hamiltonian 



1/2 



L -L 

2x^tot 



(65) 



where the superscript reminds us that the lowest allowed 
angular momentum quantum number / of the modified 
angular momentum operator L is / 1/2. 

In the notation of Ref [l^, the angular wavefunction of 
the I = 1/2, m; = 1/2 ground state of this modified rotor 
Hamiltonian is the monopole harmonic ^1/2, 1/2, 1/2(^1 0)- 
(63) To model (n^)^, we must compute the expectation value 



cos(0)) 1/2,1/2,1/2 Bund multiply this result by m 
To do this we note that 



\Yi/2,i/2,±i/2{0, <^)P = 4^(1 ± cos(0)) , (66) 
which immediately implies 



(n")t = m J dcos(6')d</>cos(6')|yi/2,i/2,i/2(6',<^)r = 

(67) 

Thus, a more general phenomenological approach that 
goes beyond sublattice-spin mean-field theory but ignores 
all non-zero wavevector modes also gives 



, m 
n = — . 
3 



(68) 



Since our QMC data show clear deviatons from this re- 
sult, we conclude that such non-zero wavevector modes 
are essential for a correct calculation of the universal 
function n^(m). 



V. DISCUSSION 

A natural question that arises from our results is 
whether the universal ground state spin texture we have 
found here can be successfully described using an effective 
field theory approach of the type used recently by Eggert 
and collaborators for studying universal aspects of the 
alternating order induced by missing spins in two dimen- 
sional S* = 1/2 antiferromagnets This approach uses 
a non-linear sigma-model description of the local anti- 
ferromagnetic order parameter, with lattice scale physics 
only entering via the values of the stiffness constant pg 
and the transverse susceptibility ^.l , and the presence of 
the vacancy captured by a local term in the action. An 
analogous treatment for our situation would need two 
things — one is a way of restricting attention to averages 




FIG. 9: Fourier transform (with antiperiodic boundary conditions assumed for convenience) of the spin- wave result for ^"{f) 
(assuming 5* = 3/2 and calculated using L — 75 for JJ2 and J J' model) along cuts passing through the antiferromagnetic 
wavevector {Tv,n). Note the nearly universal nature of the results in the neighbourhood of the antiferromagnetic wavevector, 
which in any case accounts for most of the weight of the transformed signal. 



in the S'tot = 1/2 component \G)i- of the ground state 
doublet, and the other is an understanding of the right 
boundary conditions or boundary terms in the action, so 
as to correctly reflect that fact that our finite sample has 
open boundaries. We leave this as an interesting direc- 
tion for future work, which may shed some light on the 
role of non-zero wavevector modes that were left out of 
the rotor description of the earlier section. 
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